rm(list = ls())
gc()
require(tidyverse)
# require(marginaleffects)
# require(patchwork)
require(mcp)

################################################################################
##
## Date:    2024-12-22
## Author:  james.h.bisbee@vanderbilt.edu
## Purpose: This script generates SI Figure 6.
## Inputs:  /scratch/jhb362/zilinsky_2023/data/results/VIMP_ranger/VIMP_2024_outcome-.*(CUTBACKSPEND|MON|FEEL|SPEND).*-weeks.*_DumFacts-FALSE_temp-chg12.RData
##            - Variable importance results generated by NFR_vimp_prep.R
##            - Summarized on the NYU HPC into PSRM_design_subset.RData via NFR_data_prep.R
## Outputs: ./figures/figS6.pdf
## Notes:   For this to run, JAGS version 4 is required.
##
################################################################################

# Compute details
print(paste0('Compute environment from ',Sys.Date(),' run by Bisbee'))
if(Sys.info()['sysname'] == 'Windows') {
  ram_size = system("wmic MemoryChip get Capacity", intern = TRUE)[-1]
  model_name = system("wmic cpu get name", intern = TRUE)[2] # nocov
  vendor_id = system("wmic cpu get manufacturer", intern = TRUE)[2] # nocov
  
  print(list(ram = stringr::str_squish(ram_size)[1],
             vendor_id = stringr::str_squish(vendor_id),
             model_name = stringr::str_squish(model_name),
             no_of_cores = parallel::detectCores()))
} else if(Sys.info()['sysname'] == 'Linuxs') {
  splitted <- strsplit(system("ps -C rsession -o %cpu,%mem,pid,cmd", intern = TRUE), " ")
  df <- do.call(rbind, lapply(splitted[-1], 
                              function(x) data.frame(
                                cpu = as.numeric(x[2]),
                                mem = as.numeric(x[4]),
                                pid = as.numeric(x[5]),
                                cmd = paste(x[-c(1:5)], collapse = " "))))
  df
} else {
  cat("If not on Linux or Windows, you'll have to figure out your own solution to seeing the compute environment.")
}

sessionInfo()

load('./data/VIMP_ranger/PSRM_design_subset.RData')

lookup <- vimp %>%
  count(outcome)

lookup$labs <- c("At this time, are you cutting back on how\nmuch money you spend each week, or not?",
                 "Have there been times in the past twelve months when you did not\nhave enough money to buy food that you or your family needed?",
                 "Are you feeling better about your\nfinancial situation these days, or not?",
                 "Are you feeling pretty good these days about the\namount of money you have to spend, or not?",
                 "[AGREE/DIS] You have more than enough\nmoney to do what you want to do.",
                 "In the last seven days, you have worried about money.",
                 "Did you worry yesterday that you\nspent too much money, or not?",
                 "[AGREE/DIS] You are watching your\nspending very closely.")

toplot <- vimp %>%
  mutate(relVimp = vimp / predErr) %>%
  select(vars,outcome,period,n,relVimp) %>%
  group_by(vars,outcome,period) %>%
  summarise_all(list(mean = ~mean(.,na.rm=T),
                     lb = ~quantile(.,.005,na.rm=T),
                     ub = ~quantile(.,.995,na.rm=T))) %>%
  ungroup() %>%
  left_join(lookup)

model = list(
  relVimp_mean ~ 1,
  ~ 1 + period2
)

plts <- list()
for(y in toplot %>% count(outcome) %>% pull(outcome)) {
  # stop()
  fit = mcp(model,data = toplot %>%
              mutate(period = as.Date(period)) %>%
              filter(period >= as.Date('2010-01-01'),
                     period <= as.Date('2016-01-01')) %>%
              filter(vars == 'PARTY',
                     outcome == y) %>%
              arrange(period) %>%
              mutate(period2 = as.numeric(period)) %>%
              data.frame(),
            par_x = "period2")
  
  plts[[y]] <- plot(fit) + 
    scale_x_continuous(labels = function(x) as.Date(x,origin = '1970-01-01'),
                       breaks = as.numeric(as.Date(c('2011-01-01','2013-01-01','2015-01-01')))) + 
    geom_vline(xintercept = as.numeric(as.Date('2013-01-01')),
               linetype = 'dashed')
}


pdf('./figures/figS6.pdf',width = 8,height = 8)
((plts$CUTBACKSPEND & 
    labs(x = NULL,y = NULL,
         title = "Party ID Variable Importance by Month",
         subtitle = 'Cutting back on spending') & 
    theme_bw()) + (plts$FEELBETTERFIN & 
                     labs(x = NULL,y = NULL,
                          subtitle = 'Feeling better about finances') & 
                     theme_bw())) / 
  ((plts$FEELGOODMON & 
      labs(x = NULL,y = 'Party ID Variable Importance',
           subtitle = 'Feeling good about money') & 
      theme_bw()) + (plts$MONMORETHANENOUGH & 
                       labs(x = NULL,y = NULL,
                            subtitle = 'More than enough money') & 
                       theme_bw())) / 
  ((plts$MONWORRYYEST & 
      labs(x = 'Date',y = NULL,
           subtitle = 'Worried about money') & 
      theme_bw()) + (plts$WATCHSPEND & 
                       labs(x = 'Date',y = NULL,
                            subtitle = 'Watching spending') & 
                       theme_bw()))
dev.off()

# EOF